function [bhat1,g,alpha]=inv_qr(y,D,X,Z,tau,A)
% [bhat1,g,alpha,vc1,betas]=inv_qr(y,D,X,Z,tau,A)
% Note, all inputs are required.  If there are no covariates or instruments, pass [].
%  y - dependent variable
%  D - endogenous RHS variable (NOTE: This implementation requires dim(D) = 1)
%  X - covariates
%  Z - excluded instruments (NOTE: This implementation takes the linear projection of D on {X,Z}
%      and uses that as the instrument.)
%  tau - quantile index (tau in (0,1))
%  A - grid search region.  If A is not specified the grid search region is constructed
%      with centered on the 2SQR region.  Note that A must be a column vector.
% Output.
%  bhat1 - coefficient estimates. 
%  g - norm(gamma(alpha))
%  alpha - alpha.
%      NOTE: After estimation, plot(alpha,g) graphs objective function.

 
if nargin==5,
    A=[];
end

X=[X,ones(size(y))];

x=[D,X];
z=[Z,X];

if isempty(A),              % Calculate grid search region using 2SQR if no region is supplied.
    dhat=z*pinv(z'*z)*(z'*D);
    PSI=[dhat,X];
    two_sqr=rq(PSI,y,tau);
    e1=y-x*two_sqr;
    mu1=mean(e1);
    sigma1=var(e1);
    vc1=([tau*(1-tau)]/(normpdf(norminv(tau,mu1,sqrt(sigma1)),mu1,sqrt(sigma1)))^2)*pinv(PSI'*PSI); % Add a Pseudo Inverse Here!
    std_2sqr=sqrt(diag(vc1));
    alpha=[two_sqr(1,1)-2*std_2sqr(1,1):std_2sqr(1,1)/20:two_sqr(1,1)+2*std_2sqr(1,1)]';
    clear PSI two_sqr e1 mu1 sigma1 vc1 std_2srq 
else,
    dhat=z*inv(z'*z)*(z'*D);
%    dhat=Z;
    alpha=A;
end

z = [dhat,X];           % Use dhat (projection of D on {X,Z} as instrument).

for i=1:size(alpha,1);
    disp(i);
    betas(:,i)=rq(z,y-alpha(i,1)*D,tau);
    g(i,1)=norm([betas(1:size(dhat,2),i)]);
end

I=find(g==min(g));
param1=alpha(I);
est1=betas(size(dhat,2)+1:size(betas,1),I);
hat1=[param1;est1];
bhat1 = hat1;
% PSI=[dhat,X];
                                    % Calculate naive standard errors based on normality
                                    % and homoskedasticity.
% e1=y-x*hat1;
% mu1=mean(e1);
% sigma1=var(e1);
% vc1=([tau*(1-tau)]/(normpdf(norminv(tau,mu1,sqrt(sigma1)),mu1,sqrt(sigma1)))^2)*inv(PSI'*PSI);

% std_err1=sqrt(diag(vc1));
% bhat1=[hat1,std_err1];

